Polaron versus bipolaron in conducting polymers: 
a density matrix renormalization group study 
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Competition between polaron and bipolaron in conjugated polymers with nondegenerate ground 
state is systematically studied in the extended Hubbard-Peierls model with the symmetry-breaking 
Brazovskii-Kirova term, using the density matrix renormalization group method combined with 
lattice optimization in the adiabatic approximation. We demonstrate that the relative stability 
of a bipolaron over two separated polarons sensitively depends on both on-site Hubbard U and 
nearest-neighbor repulsion V. When U is much larger than V, the bipolaron state is more stabilized 
compared with mean field calculations. 

PACS numbers: 71.10.Fd,71.10.-w,71.20.Rv,71.38.+i 
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The formation of a bipolaron, in which two electrons 
are trapped in a self-induced lattice deformation, is an 
important subject in many areas of condensed matter 
physics, such as amorphous semiconductors Q, insulat- 
ing crystals Q), high temperature superconductors H|}, 
and low-dimensional materials |^]. A simple model of 
localized linear electron-lattice coupling gives us a crite- 
rion for bipolaron formation as that the electron-lattice 
coupling is strong enough to surpass Coulomb repulsion 
between the two electrons — the situation being called 
"negative effective J7" Q. However, real situations are 
more delicate in the sense that bipolarons possess internal 
degrees of freedom to allow modification of wavefunctions 
and lattice deformations so as to reduce Coulomb repul- 
sion and yet to increase binding by lattice distortion. In 
addition, a purely electronic mechanism of negative ef- 
fective U has been proposed ||. 

In one-dimensional conjugated polymers, polarons and 
bipolarons belong to the category of solitonic nonlinear 
excitations @. Solitons can exist in polymers with de- 
generate ground states, such as in trans-polyacetylene, 
whereas polarons and bipolarons are candidates for major 
charge carriers in polymers without ground state degen- 
eracy, such as in polythiophene and poly(p-phyenylene 
vinylene) . Experimentally, the existence of bipolarons in 
doped conjugated polymers has been controversial 
The initial assignments of bipolarons relied upon the 
observation of two peaks in doping-induced absorption 
spectra |J. However, theoretical calculations [[10 11 as 
well as experimental studies of doped oligomers p |Io| 
have established that the two-peak feature actually cor- 
responds to a polaron rather than a bipolaron. On the 
other hand, the existence of bipolarons has been sug- 
gested by the observation of a single photoinduced ab- 
sorption peak in an improved sample of polythiophene 

On the theoretical side, models with neglecting 
electron-electron (e-e) interactions predict that a bipo- 
laron (BP) is always energetically more stable than a 



separated pair of polarons (2P) JIJ]. The effect of long- 
range e-e interaction was studied systematically within 
the unrestricted Hartree-Fock (UHF) approximation jl5| , 
and it was shown that the interaction significantly sup- 
presses the stability of BP. A recent study in the ex- 
tended Hubbard-Peierls model by using the density ma- 
trix renormalization group (DMRG) method Jll| has 
shown that BP is favored over 2P even for a very large 
on-site repulsion U [^7). However, off-site interactions 
were not taken into account in that study, and, in view 
of the spatial extensions of a polaron and a bipolaron, it 
is very important to take those interactions as well. 

In the present paper, we investigate this problem sys- 
tematically by using the DMRG method and by taking 
both on-site U and nearest neighbor repulsion V into 
account. The DMRG method developed originally by 
White ||l6|| for quantum spin systems is quite a powerful 
method to treat strongly correlated electron systems in 
one dimension. For the present purpose we combine the 
DMRG method with lattice optimization by using the 
Hellmann-Fcynman force equilibrium condition. We will 
demonstrate that higher order correlation effects beyond 
the mean field approximation are crucial for the forma- 
tion of BP and that the stability of BP is very sensitive 
to both U and V. 

The model Hamiltonian is the one-dimensional ex- 
tended Hubbard-Peierls Hamiltonian, 

H = -^U. l+ i[c\ s c l+l . s + h.c] 

i,s 

+U J~] n^iiii + V^(rii - l)(n l+1 - 1) 
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Here, the creation operator of a 7r electron with spin s 
at site i is denoted by ct , and the number operators are 

defined as rij lS = cj s Ci. s and n, = X) s n M- The transfer 
integral, ti i+i, is assumed in the form 



1 



t 



-iySo]t + ayi 



where So is the Brazovskii-Kirova symmetry-breaking j 
rameter to lift the degeneracy of the ground state [JT] 
Ui the change of bond length between z-th and (i + '. 
th sites, and a the electron-lattice coupling constant 
the Su-Schrieffer-Heeger type The parameter K 

the spring constant for the change of bond length, a 
the last term in eq. ([!]) is added to keep the chain leng 
constant. We take the dimensionless electron-lattice cc 
pling constant A = 2a 2 /irtK = 0.136 and a/K = 0.08' 
as in Ref . [|l5| . Although the relative stability of BP a 
2P depends on the electron-lattice coupling, we fix tl 
parameter and vary e-e interaction strengths U and 
for clarity. 

We apply the finite-system DMRG algorithm 
treat the electronic ground state. The open bounda 
condition is imposed to the system with ./V sites. The 1; 
tice is treated within the adiabatic approximation. T 
bond variable yt is determined by using the Hellmar 
Feynman force equilibrium condition [|15lJl7|] . 
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FIG. 1. Optimized bond alternation patterns of a bipo- 
laron and a pair of polarons calculated in the DMRG (solid 
line) and UHF (broken line) methods for U = IM, V/U = 0.5, 
and S = 0.05. 



where Pij+i = J2s < C L C «+M > witn < " * ' > den< 
ing the ground state expectation value. The conditi 
11^ = gives T = -{2a/ N - 1) £ m j) ra , m+1 . Thus t 
lattice configuration {yi} and the electronic ground stc 
must be determined in a self-consistent fashion. We i 
an iterative method as follows: (i) Set an initial latt: 
configuration {y^}- (fi) Calculate the electronic grou 
state by the infinite-system DMRG method for {yf- j-. 
(iii) Carry out the finite-size DMRG procedure to refine 
the ground state, (iv) Calculate new lattice configura- 
tion {y^} from eq. (||). (v) Replace {y^ } with {y^} 
and go to step (ii). The procedure is continued until the 
difference between } and } becomes negligibly 
small. In the present DMRG calculations, we truncate 
the Hilbert space by keeping m = 80 states in each block, 
which turned out to be sufficient to achieve negligible in- 
dependence of calculated results. We have also checked 
the accuracy of our DMRG calculations by comparing 
with the results of numerical exact diagonalization up to 
12 sites, with the exact solution of the Hubbard model, 
and with exact numerical results in the non-interacting 
case (U = V= 0). 

We consider the doping of two electrons in a finite chain 
to examine the formation of BP or 2P. We present results 
for the system size N — 120 and the total electron num- 
ber N c = 122. Here electrons with up and down spins are 
added to the half-filled system so as to give S z — 0. To 
obtain a 2P solution, we can equally use two up spin elec- 
trons with S z = 1 instead, because the interference be- 
tween the two well-separated polarons is negligible. This 
saves a lot of CPU time especially when 2P is less stable 



In Fig. 1 we show an example of obtained lattice dis- 
tortions of BP and 2P along with corresponding UHF 
results for comparison. In this example with U = 1.5t 
and V/U = 0.5, we see that polaronic distortions (es- 
pecially for a bipolaron) are wider and shallower in the 
DMRG than in the UHF. However, if we go to much 
smaller V/U, this tendency is reversed, i.e., the distor- 
tion becomes narrower and deeper in the DMRG. This 
implies that higher order correlation beyond the UHF 
works in opposite directions with respect to U and V. 



Figure 2 displays the energy difference AE between 
the obtained BP and 2P states as a function of U with 
constant V/U for the case So = 0.05. AE > implies 
that BP is more stable than 2P. The UHF results are 
also depicted for comparison. In the case of V/U = 0.5, 
BP is more stable than 2P for U below a critical value 
U c ~ 1.8<. The allover behavior is qualitatively similar 
to the UHF result, although U c is substantially larger. 
We could not find out a metastable BP state for U > U c 
in the DMRG calculation in contrast to the UHF result. 
In this case the iteration procedure always leads to a 2P 
state with S z = even if we start from a BP config- 
uration. This is the reason why the DMRG curve for 
V/U = 0.5 in Fig. 2 is terminated at U c above which the 
2P state is the only self-consistent solution. 
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FIG. 2. Energy difference between the bipolaron 
two-polaron states against U calculated by the DA 
method (solid lines) and the UHF approximation (br< 
lines) for V/U = 0.5 and 0.25 with <5 =0.05. 
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FIG. 3. Ground state phase diagram of the two-electron 
doped system with V/U = 0.5. The solid (broken) line indi- 
cates the phase boundary between the regions of a bipolaron 
(BP) and two separated polarons (2P) obtained by the DMRG 
method (the UHF approximation). 



On the other hand, AE for V/U = 0.25 in Fig 
is always positive and does not show any indicatio 
crossing AE — even for much larger U. That is, 
is always stable in this case. This is completely diffe 
from the UHF result, which is not very sensitive to \ 
We note that this result is consistent with Ref. 17 
the limiting case of V — 0. 



The destabilization of the BP state with increasing U 
is understandable at the UHF level as follows pq |. In 
restricting our consideration to the effects of U, a bipo- 
laron with spatial extension £ B p costs energy ~ (7/£ BP , 
while a polaron with an extension £p hardly costs en- 
ergy because of largely spin-polarized excess charges in 
the case of large U. Higher order correlation may reduce 
the probability of double occupancy, thus weakening the 
[/-dependence of AE in the DMRG method. 



We have carried out similar calculations for various So . 
By plotting the critical value U c for BP destabilization 
against <5o, we obtain the ground state phase diagram in 
the parameter space of So and U. Figure 3 is the obtained 
phase diagram for V/U = 0.5, where the solid (broken) 
line indicates the phase boundary between the regions 
of stable BP and 2P in the DMRG (UHF) results. The 
phase boundary becomes lower with raising 5q, because 
the increment of So results in the decrement of £bp an d 
£p, thereby enhancing the effect of Coulomb repulsion. 
Higher order correlation stabilizes the BP state, as men- 
tioned above. In the case of V/U = 0.25 we did not find 
out the region of 2P being more stable than BP for any 
So in the DMRG method. 



In order to clarify the role of higher order correlation 
beyond the UHF approximation for the energies of 2P 
and BP, we h^ya evajluate^l ^Jd^jc^u^tigij energies in the 
obtained DMRG wave functions. The spin fluctuation 
energy, J7£\ < ( n ^~ < n i\ >)(«a~ < n U >) >, is al- 
ways larger for 2P than in BP, destabilizing the 2P state. 
Its difference between 2P and BP grows with increasing 
U, so that this contribution is more important for smaller 
V/U. On the other hand, the charge fluctuation energy, 
v 12i, s , s > < ( n i,s- < rii.a >)(n l+1 . s >- < n i+ i <s > >) >, is 
almost the same for 2P and BP, so that the net effect to 
the relative stability is negligible. Therefore, the main 
difference between the DMRG and UHF results comes 
from spin fluctuations associated with on-site Coulomb 
repulsion U. 

In summary, we have investigated the effect of e-c 
interactions on the spatial profiles and energies of po- 
larons and bipolarons in conjugated polymers by using 
the DMRG method. We have demonstrated that the 
relative stability of a bipolaron over a pair of separated 
polarons is sensitive to both on-site U and nearest neigh- 
bor V. Spin fluctuations associated with on-site U con- 
tributes to the stabilization of a bipolaron compared with 
the results of the UHF approximation. On the other 
hand, nearest neighbor V works in the direction of dis- 
sociating a bipolaron into two polarons in the same way 
as the UHF calculation. 

The present study is limited to the extended Hubbard- 
Peierls model with short range interactions, while longer 
range interactions such as the Ohno potential should be 
taken into account to discuss situations in real materi- 
als. Judging from the roles of U and V clarified in the 
present work, we expect that the inclusion of longer range 
interactions would not cause qualitative differences but 
enhance the tendency of stabilizing polarons by V. 
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